%-------------------------------------------------------------------------------

% This file is part of code_saturne, a general-purpose CFD tool.
%
% Copyright (C) 1998-2025 EDF S.A.
%
% This program is free software; you can redistribute it and/or modify it under
% the terms of the GNU General Public License as published by the Free Software
% Foundation; either version 2 of the License, or (at your option) any later
% version.
%
% This program is distributed in the hope that it will be useful, but WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
% details.
%
% You should have received a copy of the GNU General Public License along with
% this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
% Street, Fifth Floor, Boston, MA 02110-1301, USA.

%-------------------------------------------------------------------------------

\programme{cfmsvl}
%
\vspace{1cm}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Fonction}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Pour les notations et l'algorithme dans son ensemble,
on se reportera à \fort{cfbase}.

On considère un premier pas fractionnaire au cours duquel l'énergie totale
est fixe. Seules varient la masse volumique et le flux de masse acoustique
normal aux faces (défini et calculé aux faces).

On a donc le système suivant, entre $t^n$ et $t^*$~:
\begin{equation}\label{Cfbl_Cfmsvl_eq_acoustique_cfmsvl}
\left\{\begin{array}{l}

\displaystyle\frac{\partial\rho}{\partial t}+\divs{\vect{Q}_{ac}} = 0 \\
\\
\displaystyle\frac{\partial\vect{Q}_{ac}}{\partial t}+\gradv{P} =
\rho \vect{f}\\
\\
\vect{Q}^*=\vect{Q}^n\\
\\
e^*=e^n\\

\end{array}\right.
\end{equation}

Une partie des termes sources de l'équation de la
quantité de mouvement peut être prise en compte dans cette étape
(les termes les plus importants, en prêtant attention aux sous-équilibres).

Il faut noter que si $\vect{f}$ est effectivement nul, on aura bien un
système ``acoustique'', mais que si l'on place des termes supplémentaires
dans $\vect{f}$, la dénomination est abusive (on la conservera cependant).

On obtient $\rho^* = \rho^{n+1}$ en résolvant (\ref{Cfbl_Cfmsvl_eq_acoustique_cfmsvl}),
et l'on actualise alors le flux de masse acoustique $\vect{Q}_{ac}^{n+1}$,
qui servira pour la convection (en particulier pour la convection de
l'enthalpie totale et de tous les scalaires transportés).

Suivant la valeur de \var{IGRDPP}, on actualise éventuellement la pression, en
utilisant la loi d'état :
$$
\displaystyle P^{Pred}=P(\rho^{n+1},\varepsilon^{n})
$$

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Discrétisation}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%---------------------------------
\subsection*{Discrétisation en temps}
%---------------------------------

Le système (\ref{Cfbl_Cfmsvl_eq_acoustique_cfmsvl}) discrétisé en temps donne :
\begin{equation}\label{Cfbl_Cfmsvl_eq_acoustique_discrete_cfmsvl}
\left\{\begin{array}{l}

\displaystyle\frac{\rho^{n+1}-\rho^n}{\Delta t^n}
+ \divs{\vect{Q}_{ac}^{n+1}} = 0 \\
\\
\displaystyle\frac{\vect{Q}_{ac}^{n+1}-\vect{Q}^n}{\Delta t^n}+\gradv{P^*} =
\rho^n \vect{f}^n\\
\\
Q^*=Q^n\\
\\
e^*=e^n\\

\end{array}\right.
\end{equation}

\begin{equation}\label{Cfbl_Cfmsvl_eq_forces_supplementaires_cfmsvl}
\begin{array}{llll}
\text{avec\ }&\vect{f}^n &=& \vect{0} \\
\text{ou\ }&\vect{f}^n &=& \vect{g} \\
\text{ou même\ }&\vect{f}^n &=& \vect{f}_v
 + \displaystyle\frac{1}{\rho^n}
\left( - \divs(\vect{u} \otimes \vect{Q}) + \divv(\tens{\Sigma}^v)
 + \vect{j}\wedge\vect{B} \right)^n
\end{array}
\end{equation}

Dans la pratique nous avons décidé de prendre $\vect{f}^n=\vect{g}$~:
\begin{itemize}
  \item le terme $\vect{j}\wedge\vect{B}$ n'a pas été testé,
  \item le terme $\divv(\tens{\Sigma}^v)$ était négligeable sur les tests
        réalisés,
  \item le terme $\divs(\vect{u} \otimes \vect{Q})$ a paru déstabiliser les
        calculs (mais au moins une partie des tests a été réalisée
        avec une erreur de programmation et il faudrait donc les reprendre).
\end{itemize}
\bigskip

Le terme $\vect{Q}^n$ dans la 2\textsuperscript{ème} équation
de (\ref{Cfbl_Cfmsvl_eq_acoustique_discrete_cfmsvl}) est le vecteur ``quantité de mouvement''
qui provient de l'étape de résolution de la quantité de mouvement du pas
de temps précédent, $\vect{Q}^n = \rho^n \vect{u}^n$.
On pourrait théoriquement utiliser un vecteur quantité de mouvement issu
de l'étape acoustique du pas de temps précédent, mais il ne constitue
qu'un ``prédicteur'' plus ou moins satisfaisant (il n'a pas ``vu'' les termes
sources qui ne sont pas dans  $\vect{f}^n$) et cette solution
n'a pas été testée.

\bigskip
On écrit alors la pression sous la forme~:
\begin{equation}
\gradv{P}=c^2\,\gradv{\rho}+\beta\,\gradv{s}
\end{equation}

avec $c^2 = \left.\displaystyle\frac{\partial P}{\partial \rho}\right|_s$
et $\beta = \left.\displaystyle\frac{\partial P}{\partial s}\right|_\rho$
tabulés ou analytiques à partir de la loi d'état.

On discrétise l'expression précédente en~:
\begin{equation}
\gradv{P^*}=(c^2)^n\gradv(\rho^{n+1})+\beta^n\gradv(s^n)
\end{equation}

On obtient alors une équation
portant sur $\rho^{n+1}$ en substituant l'expression de $\vect{Q}_{ac}^{n+1}$
issue de la 2\textsuperscript{ème} équation
de~(\ref{Cfbl_Cfmsvl_eq_acoustique_discrete_cfmsvl})
dans la 1\textsuperscript{ère} équation
de~(\ref{Cfbl_Cfmsvl_eq_acoustique_discrete_cfmsvl})~:
\begin{equation}\label{Cfbl_Cfmsvl_eq_densite_cfmsvl}
\displaystyle\frac{\rho^{n+1}-\rho^n}{\Delta t^n}
+\divs(\vect{w}^n \rho^n)
-\divs\left(\Delta t^n (c^2)^n \gradv(\rho^{n+1})\right) = 0
\end{equation}

où~:
\begin{equation}
\begin{array}{lll}
\vect{w}^n&=&  \vect{u}^n + \Delta t^n
\displaystyle\left(\vect{f}^n-\frac{\beta^n}{\rho^n}\gradv(s^n)\right)
\end{array}
\end{equation}

Formulation alternative (programmée mais non testée)
avec le terme de convection implicite~:
\begin{equation}\label{Cfbl_Cfmsvl_eq_densite_bis_cfmsvl}
\displaystyle\frac{\rho^{n+1}-\rho^n}{\Delta t^n}
+\divs(\vect{w}^n \rho^{n+1})
-\divs\left(\Delta t^n (c^2)^n \gradv(\rho^{n+1})\right) = 0
\end{equation}


%---------------------------------
\subsection*{Discrétisation en espace}
%---------------------------------


%.................................
\subsubsection*{Introduction}
%.................................

On intègre l'équation précédente ( (\ref{Cfbl_Cfmsvl_eq_densite_cfmsvl})
ou (\ref{Cfbl_Cfmsvl_eq_densite_bis_cfmsvl}) ) sur la cellule $i$ de volume $\Omega_i$.
On transforme les intégrales de volume en intégrales surfaciques
et l'on discrétise ces intégrales. Pour simplifier l'exposé, on se
place sur une cellule $i$ dont aucune face n'est sur le bord du domaine.

On obtient alors l'équation discrète
suivante\footnote{L'exposant $^{n+\frac{1}{2}}$ signifie que le terme
peut être implicite ou explicite. En pratique on a choisi
$\rho^{n+\frac{1}{2}} = \rho^{n}$.}~:
\begin{equation}\label{Cfbl_Cfmsvl_eq_densite_discrete_cfmsvl}
\Omega_i \displaystyle\frac{\rho_i^{n+1}-\rho_i^n}{\Delta t^n}
+\sum\limits_{j\in Vois(i)}(\rho^{n+\frac{1}{2}} \vect{w}^n)_{ij} \cdot \vect{S}_{ij}
-\sum\limits_{j\in Vois(i)} \left(\Delta t^n (c^2)^n
\gradv(\rho^{n+1})\right)_{ij} \cdot \vect{S}_{ij}
= 0
\end{equation}

%.................................
\subsubsection*{Discrétisation de la partie ``convective''}
%.................................

La valeur à la face s'écrit~:
\begin{equation}
(\rho^{n+\frac{1}{2}} \vect{w}^n)_{ij} \cdot \vect{S}_{ij}
= \rho^{n+\frac{1}{2}}_{ij} \vect{w}^n_{ij} \cdot \vect{S}_{ij}
\end{equation}
avec, pour $\vect{w}^n_{ij}$,
une simple interpolation linéaire~:
\begin{equation}
\vect{w}^n_{ij}
= \alpha_{ij} \vect{w}^n_i + (1-\alpha_{ij}) \vect{w}^n_j
\end{equation}
et un décentrement sur la valeur de $\rho^{n+\frac{1}{2}}$ aux faces~:
\begin{equation}
\begin{array}{lllll}
\displaystyle\rho_{ij}^{n+\frac{1}{2}} &=& \rho_{I'}^{n+\frac{1}{2}}
                &\text{si\ }& \vect{w}^n_{ij} \cdot \vect{S}_{ij} \geqslant 0 \\
                         &=& \rho_{J'}^{n+\frac{1}{2}}
                &\text{si\ }& \vect{w}^n_{ij} \cdot \vect{S}_{ij} < 0 \\
\end{array}
\end{equation}
que l'on peut noter~:
\begin{equation}
\displaystyle\rho_{ij}^{n+\frac{1}{2}}
 = \beta_{ij}\rho_{I'}^{n+\frac{1}{2}} + (1-\beta_{ij})\rho_{J'}^{n+\frac{1}{2}}
\end{equation}
avec
\begin{equation}
\left\{\begin{array}{lll}
\beta_{ij} = 1 & \text{si\ } & \vect{w}^n_{ij} \cdot \vect{S}_{ij} \geqslant 0 \\
\beta_{ij} = 0 & \text{si\ } & \vect{w}^n_{ij} \cdot \vect{S}_{ij} < 0 \\
\end{array}\right.
\end{equation}

%.................................
\subsubsection*{Discrétisation de la partie ``diffusive''}
%.................................

La valeur à la face s'écrit~:
\begin{equation}
\left(\Delta t^n (c^2)^n \gradv(\rho^{n+1})\right)_{ij}\cdot \vect{S}_{ij}
= \Delta t^n (c^2)^n_{ij}
\displaystyle \left( \frac{\partial \rho}{\partial n} \right)^{n+1}_{ij}S_{ij}
\end{equation}
avec, pour assurer la continuité du flux normal à l'interface,
une interpolation harmonique de $(c^2)^n$~:
\begin{equation}\label{Cfbl_Cfmsvl_eq_harmonique_cfmsvl}
\displaystyle(c^2)_{ij}^n
= \frac{(c^2)_{i}^n (c^2)_{j}^n}
{\alpha_{ij}(c^2)_{i}^n+(1-\alpha_{ij})(c^2)_{j}^n}
\end{equation}
et un schéma centré pour le gradient normal aux faces~:
\begin{equation}
\displaystyle \left( \frac{\partial \rho}{\partial n} \right)^{n+1}_{ij}
= \displaystyle\frac{\rho_{J'}^{n+1}-\rho_{I'}^{n+1}}{\overline{I'J'}}
\end{equation}

%.................................
\subsubsection*{Système final}
%.................................

On obtient maintenant le système final, portant sur
$(\rho_i^{n+1})_{i=1 \ldots N}$~:
\begin{equation}\label{Cfbl_Cfmsvl_eq_densite_finale_cfmsvl}
\displaystyle\frac{\Omega_i}{\Delta t^n} (\rho_i^{n+1}-\rho_i^n)
+\sum\limits_{j\in Vois(i)}\rho_{ij}^{n+\frac{1}{2}}
\vect{w}_{ij}^n \cdot \vect{S}_{ij}
-\sum\limits_{j\in Vois(i)} \Delta t^n (c^2)_{ij}^n
\displaystyle\frac{\rho_{J'}^{n+1}-\rho_{I'}^{n+1}}{\overline{I'J'}}\ S_{ij}
= 0
\end{equation}



%.................................
\subsubsection*{Remarque~: interpolation aux faces pour le terme de diffusion}
%.................................

Le choix de la forme de la moyenne pour le cofacteur du flux
normal n'est pas sans conséquence sur la vitesse de convergence, surtout
lorsque l'on est en présence de fortes inhomogénéités.

On utilise une interpolation harmonique pour $c^2$
afin de conserver la continuité du flux diffusif normal
$\Delta t (c^2) \displaystyle\frac{\partial \rho}{\partial n}$
à l'interface $ij$. En effet, on suppose que le flux est dérivable à
l'interface. Il doit donc y être continu.\\
%
Écrivons la continuité du flux normal à l'interface,
avec la discrétisation
suivante\footnote{On ne reconstruit pas les valeurs de $\Delta\,t\,c^2$
aux points $I'$ et
$J'$.}~:
\begin{equation}
\left(\Delta t (c^2)\displaystyle\frac{\partial \rho}{\partial n}\right)_{ij}
= \Delta t (c^2)_i  \displaystyle\frac{\rho_{ij} - \rho_{I'}   }{\overline{I'F}}
=  \Delta t (c^2)_j  \displaystyle\frac{\rho_{J'}    - \rho_{ij}}{\overline{FJ'}}
\end{equation}
En égalant les flux à gauche et à droite de l'interface, on obtient
\begin{equation}
\rho_{ij} = \displaystyle\frac{\overline{I'F}\,(c^2)_j\rho_{J'} + \overline{FJ'}\,(c^2)_i\rho_{I'}}
{\overline{I'F}\,(c^2)_j + \overline{FJ'}\,(c^2)_i}
\end{equation}
On introduit cette formulation dans la définition du flux (par exemple, du
flux à gauche)~:
\begin{equation}
\left(\Delta t (c^2)\displaystyle\frac{\partial \rho}{\partial n}\right)_{ij}
= \Delta t (c^2)_i  \displaystyle\frac{\rho_{ij} - \rho_{I'}   }{\overline{I'F}}
\end{equation}
et on utilise la définition de $(c^2)_{ij}$ en fonction de ce même flux
\begin{equation}
\left(\Delta t (c^2)\displaystyle\frac{\partial \rho}{\partial n}\right)_{ij}
 \stackrel{\text{déf}}{=}
 \Delta t (c^2)_{ij} \displaystyle\frac{\rho_{J'}    - \rho_{I'}   }{\overline{I'J'}}
\end{equation}
pour obtenir la valeur de $(c^2)_{ij}$ correspondant à l'équation (\ref{Cfbl_Cfmsvl_eq_harmonique_cfmsvl})~:
\begin{equation}
(c^2)_{ij} = \displaystyle\frac{\overline{I'J'}\,(c^2)_i(c^2)_j}{\overline{FJ'}\,(c^2)_i + \overline{I'F}\,(c^2)_j}
\end{equation}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Mise en \oe uvre}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Le système (\ref{Cfbl_Cfmsvl_eq_densite_finale_cfmsvl}) est résolu par une méthode
d'incrément et résidu en utilisant
une méthode de Jacobi pour inverser le système si le terme convectif
est implicite et en utilisant une méthode de gradient conjugué
si le terme convectif est explicite (qui est le cas par défaut).

Attention, les valeurs du flux de masse $\rho\,\vect{w}\cdot\vect{S}$ et
de la viscosité $\Delta\,t\,c^2\frac{S}{d}$ aux faces de
bord, qui sont calculées dans \fort{cfmsfl} et \fort{cfmsvs} respectivement,
sont modifiées immédiatement après l'appel à ces sous-programmes.
En effet, il est indispensable que la contribution de bord de
$\left(\rho\,\vect{w}-\Delta\,t\,(c^2)\,\gradv\,\rho\right)\cdot\vect{S}$
représente exactement $\vect{Q}_{ac}\cdot\vect{S}$.
Pour cela,
\begin{itemize}
\item immédiatement après l'appel à
\fort{cfmsfl}, on remplace la contribution de bord de
$\rho\,\vect{w}\cdot\vect{S}$
par le flux de masse exact, $\vect{Q}_{ac}\cdot\vect{S}$,
déterminé à partir des conditions aux limites,
\item puis, immédiatement après l'appel à
\fort{cfmsvs}, on annule la viscosité au bord $\Delta\,t\,(c^2)$ pour
éliminer la contribution de $-\Delta\,t\,(c^2)\,(\gradv\,\rho)\cdot\vect{S}$
(l'annulation de la viscosité n'est pas problématique pour la matrice,
puisqu'elle porte sur des incréments).
\end{itemize}

\bigskip

Une fois qu'on a obtenu $\rho^{n+1}$,
on peut actualiser le flux de masse acoustique
aux faces $(\vect{Q}_{ac}^{n+1})_{ij} \cdot \vect{S}_{ij}$,
qui servira pour la convection des autres variables~:
\begin{equation}\label{Cfbl_Cfmsvl_eq_flux_masse_acoustique_cfmsvl}
\displaystyle(\vect{Q}_{ac}^{n+1})_{ij}\cdot\vect{S}_{ij}=
-\left(\Delta t^n (c^2)^n \gradv(\rho^{n+1})\right)_{ij}\cdot\vect{S}_{ij}
+\left(\rho^{n+\frac{1}{2}} \vect{w}^n\right)_{ij}\cdot\vect{S}_{ij}\\
\end{equation}
Ce calcul de flux est réalisé par \fort{cfbsc3}.
Si l'on a choisi l'algorithme standard, équation~(\ref{Cfbl_Cfmsvl_eq_densite_cfmsvl}),
on complète le flux dans \fort{cfmsvl} immédiatement après l'appel
à \fort{cfbsc3}.
En effet, dans ce cas,
la convection est explicite ($\rho^{n+\frac{1}{2}}=\rho^{n}$,
obtenu en imposant \var{ICONV(ISCA(IRHO))=0})
et le sous-programme \fort{cfbsc3},
qui calcule le flux de masse aux faces,
ne prend pas en compte la contribution du terme
$\rho^{n+\frac{1}{2}}\,\vect{w}^n\cdot\vect{S}$. On ajoute donc cette
contribution dans \fort{cfmsvl}, après l'appel à \fort{cfbsc3}.
Au bord, en particulier, c'est bien le flux de masse calculé à partir
des conditions aux limites que l'on obtient.

On actualise la pression à la fin de l'étape, en utilisant la loi d'état~:
\begin{equation}
\displaystyle P_i^{pred}=P(\rho_i^{n+1},\varepsilon_i^{n})
\end{equation}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Points à traiter}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Le calcul du flux de masse au  bord n'est pas entièrement satisfaisant
si la convection est traitée de manière implicite
(algorithme non standard, non testé,
associé à l'équation~(\ref{Cfbl_Cfmsvl_eq_densite_bis_cfmsvl}),
correspondant au choix $\rho^{n+\frac{1}{2}}=\rho^{n+1}$ et
obtenu en imposant \var{ICONV(ISCA(IRHO))=1}).
En effet, après \fort{cfmsfl}, il faut déterminer la vitesse de
convection $\vect{w}^n$ pour qu'apparaisse
$\rho^{n+1} \vect{w}^n\cdot\vect{n}$
au cours de la résolution par \fort{cs\_equation\_iterative\_solve}. De ce fait, on doit déduire
une valeur de $\vect{w}^n$ à partir de la valeur
du flux de masse. Au bord, en particulier, il faut
donc diviser le flux de masse
issu des conditions aux limites par la valeur de bord de $\rho^{n+1}$.
Or, lorsque des conditions de Neumann sont appliquées à la
masse volumique,
la valeur de $\rho^{n+1}$ au bord n'est pas connue avant la résolution du
système. On utilise donc, au lieu de la valeur de bord inconnue de
$\rho^{n+1}$ la valeur de bord prise au pas de temps
précédent $\rho^{n}$. Cette approximation est susceptible
d'affecter la valeur du flux de masse au bord.
